
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

from sklearn.preprocessing import StandardScaler, LabelEncoder
from sklearn.cluster import KMeans
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import classification_report, confusion_matrix, accuracy_score, recall_score, precision_score, f1_score

import xgboost as xgb
from imblearn.over_sampling import SMOTE
from sklearn.model_selection import GridSearchCV
from imblearn.pipeline import Pipeline
from sklearn.feature_selection import SelectKBest

from scipy import stats
from statsmodels.stats.outliers_influence import variance_inflation_factor
from kneed import KneeLocator
from joblib import Parallel, delayed

import plotly.express as px
import plotly.graph_objs as go
import plotly.io as pio
pio.renderers.default = "plotly_mimetype+notebook"


import warnings
warnings.filterwarnings("ignore")


# Preprocess data based on the provided data structure
def preprocess_data(df):
    # Handle missing values
    print(list(df.columns))
    df.replace('Unknown',0,inplace=True)
    df['HandsetPrice'] = df['HandsetPrice'].astype('float')
    df.fillna(df.mean(), inplace=True)
    df.fillna(df.mode().iloc[0], inplace=True)
    
    # Drop unnecessary columns (if any)
    # df = df.drop(['ColumnToDrop'], axis=1)
    
    return df

def remove_constant_features(df, threshold=0):
    # Identify constant features
    constant_features = [col for col in df.columns if df[col].nunique() <= threshold]
    print(f"Removing constant features: {constant_features}")
    
    # Drop constant features
    df = df.drop(columns=constant_features)
    
    return df

def load_and_preprocess_data(file_path):
    # Load the data
    df = pd.read_csv(file_path)
    
    # Perform preprocessing and feature extraction
    df = preprocess_data(df)

    # Remove constant features
    df = remove_constant_features(df)
    
    # Label encode the binary categorical columns
    binary_categorical_columns = ['ChildrenInHH', 'HandsetRefurbished', 'HandsetWebCapable', 'TruckOwner', 'RVOwner',
                                  'Homeownership', 'BuysViaMailOrder', 'RespondsToMailOffers', 'OptOutMailings',
                                  'NonUSTravel', 'OwnsComputer', 'HasCreditCard', 'NewCellphoneUser',
                                  'NotNewCellphoneUser', 'OwnsMotorcycle', 'MadeCallToRetentionTeam']
    for col in binary_categorical_columns:
        le = LabelEncoder()
        df[col] = le.fit_transform(df[col].astype(str))

    # Evaluate correlation and multicollinearity before one-hot encoding
    df = evaluate_correlation_and_multicollinearity(df)

    # One-hot encode the remaining categorical columns
    categorical_columns = ['ServiceArea', 'CreditRating', 'PrizmCode', 'Occupation', 'MaritalStatus']
    df = pd.get_dummies(df, columns=categorical_columns, drop_first=True)  # Add drop_first=True

    return df


def evaluate_correlation_and_multicollinearity(df, correlation_threshold=0.8, vif_threshold=5, cramers_v_threshold=0.7):
    # Step 1: Calculate the correlation matrix for numerical variables
    numerical_cols = df.select_dtypes(include=['int64', 'float64']).columns
    #numerical_cols = numerical_cols.drop('Churn')  # Exclude 'Churn' from numerical_cols
    correlation_matrix = df[numerical_cols].corr()

    # Step 2: Use the Chi-squared test of independence for categorical variables
    categorical_cols = df.select_dtypes(include=['object']).columns
    categorical_cols = categorical_cols.drop('Churn')
    chi2_matrix = pd.DataFrame(index=categorical_cols, columns=categorical_cols)
    for col1 in categorical_cols:
        for col2 in categorical_cols:
            contingency_table = pd.crosstab(df[col1], df[col2])
            _, p_value, _, _ = stats.chi2_contingency(contingency_table)
            chi2_matrix.loc[col1, col2] = p_value

    # Step 3: Calculate the VIF to detect multicollinearity among numerical variables
    numerical_df = df[numerical_cols]#.drop(columns='Churn')
    vif = pd.DataFrame()
    vif['Feature'] = numerical_df.columns
    vif['VIF'] = [variance_inflation_factor(numerical_df.values, i) for i in range(numerical_df.shape[1])]

    # Step 4: Use Cramér's V for categorical variables
    def cramers_v(confusion_matrix):
        chi2 = stats.chi2_contingency(confusion_matrix)[0]
        n = confusion_matrix.sum().sum()
        phi2 = chi2 / n
        r, k = confusion_matrix.shape
        phi2corr = max(0, phi2 - ((k - 1) * (r - 1)) / (n - 1))
        rcorr = r - ((r - 1) ** 2) / (n - 1)
        kcorr = k - ((k - 1) ** 2) / (n - 1)
        return np.sqrt(phi2corr / min((kcorr - 1), (rcorr - 1)))

    cramers_v_matrix = pd.DataFrame(index=categorical_cols, columns=categorical_cols)
    for col1 in categorical_cols:
        for col2 in categorical_cols:
            confusion_matrix = pd.crosstab(df[col1], df[col2])
            cramers_v_value = cramers_v(confusion_matrix)
            cramers_v_matrix.loc[col1, col2] = cramers_v_value

    # Remove highly correlated numerical features
    correlated_features = set()
    for i in range(len(correlation_matrix.columns)):
        for j in range(i):
            if abs(correlation_matrix.iloc[i, j]) > correlation_threshold:
                colname = correlation_matrix.columns[i]
                correlated_features.add(colname)

    # Remove numerical features with high VIF
    high_vif_features = vif[vif['VIF'] > vif_threshold]['Feature']
    correlated_features.update
    high_vif_features.values.tolist()

    # Remove categorical features with low p-values (high association) and high Cramér's V values
    for col1 in categorical_cols:
        for col2 in categorical_cols:
            if (col1 != col2) and (chi2_matrix.loc[col1, col2] < (1 - correlation_threshold)) and (cramers_v_matrix.loc[col1, col2] > cramers_v_threshold):
                correlated_features.add(col1)

    # Drop correlated features
    df = df.drop(columns=list(correlated_features))

    return df


# Replace with the path to your customer data
file_path = "./data/cell2celltrain.csv"
df = load_and_preprocess_data(file_path)

['CustomerID', 'Churn', 'MonthlyRevenue', 'MonthlyMinutes', 'TotalRecurringCharge', 'DirectorAssistedCalls', 'OverageMinutes', 'RoamingCalls', 'PercChangeMinutes', 'PercChangeRevenues', 'DroppedCalls', 'BlockedCalls', 'UnansweredCalls', 'CustomerCareCalls', 'ThreewayCalls', 'ReceivedCalls', 'OutboundCalls', 'InboundCalls', 'PeakCallsInOut', 'OffPeakCallsInOut', 'DroppedBlockedCalls', 'CallForwardingCalls', 'CallWaitingCalls', 'MonthsInService', 'UniqueSubs', 'ActiveSubs', 'ServiceArea', 'Handsets', 'HandsetModels', 'CurrentEquipmentDays', 'AgeHH1', 'AgeHH2', 'ChildrenInHH', 'HandsetRefurbished', 'HandsetWebCapable', 'TruckOwner', 'RVOwner', 'Homeownership', 'BuysViaMailOrder', 'RespondsToMailOffers', 'OptOutMailings', 'NonUSTravel', 'OwnsComputer', 'HasCreditCard', 'RetentionCalls', 'RetentionOffersAccepted', 'NewCellphoneUser', 'NotNewCellphoneUser', 'ReferralsMadeBySubscriber', 'IncomeGroup', 'OwnsMotorcycle', 'AdjustmentsToCreditRating', 'HandsetPrice', 'MadeCallToRetentionTeam', 'CreditRating', 'PrizmCode', 'Occupation', 'MaritalStatus']
Removing constant features: []


with pd.option_context('display.max_rows', None, 'display.max_columns', None):
    display(df.head())
    display(df.info( show_counts=True))

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 51047 entries, 0 to 51046
Columns: 816 entries, CustomerID to Segment
dtypes: float64(24), int32(18), int64(9), object(1), uint8(764)
memory usage: 53.9+ MB

None


def kmeans_inertia(k, scaled_data):
    kmeans = KMeans(n_clusters=k, random_state=42)
    kmeans.fit(scaled_data)
    return kmeans.inertia_

# Perform customer segmentation
def customer_segmentation(df, max_k=10):
    X = df.drop(['CustomerID','Churn'], axis=1)
    #print(list(X.columns))
    X = StandardScaler().fit_transform(X)
    
    # Find the optimal number of clusters (K) using the Elbow method
    sum_of_squared_distances = []
    k_values = range(1, max_k+1)

    sum_of_squared_distances= (Parallel(n_jobs=-1)(delayed(kmeans_inertia)(k, X) for k in k_values))

    # Identify the elbow point programmatically using the KneeLocator
    kneedle = KneeLocator(k_values, sum_of_squared_distances, curve='convex', direction='decreasing')
    optimal_k = kneedle.knee
    print(f"Optimal K value: {optimal_k}")
    
    model = KMeans(n_clusters=optimal_k, random_state=42)
    labels = model.fit_predict(X)
    
    # Calculate the mean value of each feature for each cluster
    df['Cluster'] = labels
    feature_means_by_cluster = df.drop(['CustomerID','Churn'], axis=1).groupby('Cluster').mean()

     # Rank the features by their importance within each cluster
    top_features_by_cluster = feature_means_by_cluster.apply(lambda s: s.abs().nlargest(5).index.tolist(), axis=1)

    # Create a grouped bar chart of the top 5 features for each cluster with unique colors
    fig = go.Figure()

    colors = px.colors.qualitative.Plotly[:optimal_k]
    x_labels = []
    for cluster, top_features in top_features_by_cluster.iteritems():
        cluster_data = feature_means_by_cluster.loc[cluster, top_features]
        cluster_data_percentage = cluster_data / cluster_data.sum() * 100
        for feature, percentage in zip(top_features, cluster_data_percentage):
            fig.add_trace(go.Bar(
                x=[f"Cluster {cluster} - {feature}"],
                y=[percentage],
                name=f"Cluster {cluster} - {feature}",
                marker_color=colors[cluster]
            ))
            x_labels.append(f"Cluster {cluster} - {feature}")

    fig.update_layout(
        barmode='group',
        xaxis_title="Clusters and Top Features",
        yaxis_title="Percentage",
        title=f"Top 5 Features for Each Cluster",
        xaxis={'tickmode': 'array', 'tickvals': x_labels, 'ticktext': x_labels, 'tickangle': -45},
        showlegend=False
    )

    fig.show()

    return labels


%%time

# Run customer segmentation
labels = customer_segmentation(df)
df['Segment'] = labels

Optimal K value: 4

CPU times: total: 39.9 s
Wall time: 31.4 s


# Train churn prediction model
def train_churn_prediction_model(df):
    X = df.drop(['CustomerID', 'Churn'], axis=1)
    y = df['Churn'].map({'Yes': 1, 'No': 0})
    
    X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)
    
    # Resample the training data using SMOTE
    smote = SMOTE(random_state=42)
    X_train_resampled, y_train_resampled = smote.fit_resample(X_train, y_train)
    
    # Create a pipeline for preprocessing, feature selection, and model training
    pipeline = Pipeline([
        ('scaler', StandardScaler()),
        ('selector', SelectKBest(k=10)),
        ('classifier', xgb.XGBClassifier(use_label_encoder=False, random_state=42, eval_metric='logloss'))
    ])
    
    # Define a parameter grid for hyperparameter tuning
    param_grid = {
        'classifier__n_estimators': [100, 200, 500],
        'classifier__learning_rate': [0.01, 0.1, 0.2],
        'classifier__max_depth': [3, 6, 9],
        'classifier__subsample': [0.5, 0.8, 1],
        'classifier__colsample_bytree': [0.5, 0.8, 1],
    }
    
    # Perform a grid search using cross-validation
    grid_search = GridSearchCV(pipeline, param_grid, scoring='recall', cv=5, n_jobs=-1, verbose=1)
    grid_search.fit(X_train_resampled, y_train_resampled)
    
    # Train the model with the best hyperparameters
    best_model = grid_search.best_estimator_
    best_model.fit(X_train_resampled, y_train_resampled)
    
    # Make predictions on the test set
    y_pred = best_model.predict(X_test)

    # create the confusion matrix
    cm = confusion_matrix(actual_labels, predicted_labels)

    # create the heatmap
    heatmap = go.Heatmap(z=cm, x=['0', '1'], y=['0', '1'], colorscale='Blues')

    # create the layout
    layout = go.Layout(title='Confusion Metrix')

    # create the figure
    fig = go.Figure(data=[heatmap], layout=layout)

    # show the figure
    fig.show()

    # Print the classification report
    print(classification_report(y_test, y_pred))
    print("Best parameters found: ", grid_search.best_params_)
    print("Accuracy:", accuracy_score(y_test, y_pred))
    print("F1 Score:", f1_score(y_test, y_pred, average='weighted'))
    print("Precision:", precision_score(y_test, y_pred, average='weighted'))
    print("Recall:", recall_score(y_test, y_pred, average='weighted'))
    
    return best_model


%%time

model = train_churn_prediction_model(df)

Fitting 5 folds for each of 243 candidates, totalling 1215 fits
              precision    recall  f1-score   support

           0       0.72      0.71      0.72      7308
           1       0.30      0.32      0.31      2902

    accuracy                           0.60     10210
   macro avg       0.51      0.51      0.51     10210
weighted avg       0.60      0.60      0.60     10210

Best parameters found:  {'classifier__colsample_bytree': 0.8, 'classifier__learning_rate': 0.2, 'classifier__max_depth': 9, 'classifier__n_estimators': 500, 'classifier__subsample': 1}
Accuracy: 0.597845249755142
F1 Score: 0.6011637045432573
Precision: 0.6047614016367154
Recall: 0.597845249755142
CPU times: total: 14min 32s
Wall time: 19min 14s

	CustomerID	Churn	MonthlyRevenue	MonthlyMinutes	TotalRecurringCharge	DirectorAssistedCalls	PercChangeMinutes	PercChangeRevenues	DroppedCalls	BlockedCalls	UnansweredCalls	CustomerCareCalls	ThreewayCalls	OutboundCalls	InboundCalls	PeakCallsInOut	OffPeakCallsInOut	CallWaitingCalls	MonthsInService	UniqueSubs	ActiveSubs	Handsets	CurrentEquipmentDays	AgeHH1	AgeHH2	ChildrenInHH	HandsetWebCapable	Homeownership	BuysViaMailOrder	RespondsToMailOffers	OwnsComputer	HasCreditCard	RetentionCalls	NewCellphoneUser	NotNewCellphoneUser	IncomeGroup	AdjustmentsToCreditRating	HandsetPrice	MadeCallToRetentionTeam	ServiceArea_MILMIL414	ServiceArea_OKCTUL918	ServiceArea_PITHOM412	ServiceArea_SEAPOR503	CreditRating_3-Good	CreditRating_4-Medium	PrizmCode_Suburban	PrizmCode_Town	Occupation_Crafts	Occupation_Other	Occupation_Professional	MaritalStatus_No	MaritalStatus_Yes	Cluster	Segment
0	3000002	Yes	24.00	219.0	22.0	0.25	-157.0	-19.0	0.7	0.7	6.3	0.0	0.0	0.0	0.0	58.0	24.0	0.3	61	2	1	2.0	361.0	62.0	0.0	0	1	1	1	1	1	1	1	0	0	4	0	30.0	1	0	0	0	1	0	0	1	0	0	0	1	1	0	3	3
1	3000010	Yes	16.99	10.0	17.0	0.00	-4.0	0.0	0.3	0.0	2.7	0.0	0.0	0.0	0.0	5.0	1.0	0.0	58	1	1	2.0	1504.0	40.0	42.0	1	0	1	1	1	1	1	0	1	0	5	0	30.0	0	0	0	1	0	0	1	1	0	0	0	1	0	1	3	3
2	3000014	No	38.00	8.0	38.0	0.00	-2.0	0.0	0.0	0.0	0.0	0.0	0.0	0.3	0.0	1.3	3.7	0.0	60	1	1	1.0	1812.0	26.0	26.0	1	0	0	0	0	0	1	0	1	0	6	0	0.0	0	1	0	0	0	1	0	0	1	1	0	0	0	1	3	3
3	3000022	No	82.28	1312.0	75.0	1.24	157.0	8.1	52.0	7.7	76.0	4.3	1.3	370.3	147.0	555.7	303.7	22.7	59	2	2	9.0	458.0	30.0	0.0	0	1	1	1	1	0	1	0	1	0	6	0	10.0	0	0	0	1	0	0	1	0	0	0	1	0	1	0	1	1
4	3000026	Yes	17.14	0.0	17.0	0.00	0.0	-0.2	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	53	2	2	4.0	852.0	46.0	54.0	0	0	1	1	1	1	1	0	0	1	9	1	10.0	0	0	1	0	0	0	0	0	0	0	0	1	0	1	3	3

Customer Churn Prediction Using Machine Learning¶

I am using a sample dataset of a telecommunications company, which contains customer information, including demographic data, usage, and billing information. The target variable is Churn, which indicates if a customer has left the company.¶

Data Preprocessing¶

In this step, we load the dataset, clean the data, and preprocess it for machine learning. We remove unnecessary columns, handle missing values, and convert categorical features using one hot encoder.¶

Correlation and Multicollinearity¶

Customer Segmentation¶

Model training and evaluation¶

Conclusion¶